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I. Introduction 

The objective of this report is to review the various methods that have been 
studied in the past to allow probabilistic analysis of dynamic response for 
systems with random parameters. In general, the mechanical parameters (i.e., 
spring, damping, Joint parameters, dead zones, etc.) may not be known exactly. 
If, for example, the variations about the nominal values are very small, then 
the dynamic response would be adequately obtained deterministically. 
However, for space structures which require precise pointing, it appears that the 
variations or uncertainties about the nominal values of the structural details and 
of the environmental conditions may be too large to be considered as 
negligible. 

Thus, these uncertainties must be accounted for on some rational basis 
which we shall assume to be defined in terms of probability distributions about 
their nominal values. The quantities of concern for describing the response of 
the structure includes displacements and velocities, as well as the distributions 
of natural frequencies. The exact statistical characterization of the response 
would yield Joint probability distributions for the response variables. Since the 
random quantities will appear as coefficients, determining the exact 
distributions will be difficult at best. Thus, certain approximations will have to 
be made. There are a number of techniques that we shall discuss that are 
available even in the non-linear case. 

In the general case, the n-mass linear structural system possesses the 
dynamical description through the differential equation 

My + Cy + Ky = f(t) (l-l) 

where M, C, K are the mass, damping and stiffness matrices, and f(t) denotes 

the external excitation on the structure. 
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We shall define the vector equations (1.1) through the vector 


x 

which represents the system (l.l) as 


y 

y 


where 


x = Ax + Bf ( t) , 



( 1 . 2 ) 


(1.3) 


A = 


0 

M _1 K 


I 

M _1 C 


B = 


0 

M 


- i 


(1.4) 


The solution of (1.3) for x o (0)=x o Is 


t 

X (t) = e~ At x 0 + / e" A(t ^ r) Bf(r)dr . (1-5) 

0 

We assume that the matrix B, determined by the mass constants, is 
deterministic. Thus, the random quantities appear in matrix A. We shall 
denote the random variable in A, as X ls . . . , X m . Further, the external 
excitations f(t) may or may not be random. 


The most important topic for engineering systems is how uncertain 
parameter values influence the accuracy of system response prediction. It often 
suffices to know how these uncertainities influence the accuracy in estimating 
the values of the national frequencies and their corresponding normal modes of 
motion in a conservative system (C= 0). 


Since linear system response prediction depends upon frequency response or 
impulsive admittance, our interest will center on natural frequencies, normal 
modes, frequency response, as well as Impulse response. The methods of 
techniques that we will describe in order to pursue the various subjects are: (1) 
Liouville’s equation; (2) perturbation methods; (3) mean square approximate 
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systems; and (4) non-linear systems, with approximation by linear systems. 
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II. Liouville Equation 

In this section we derive the method based upon the Liouville equation for 
the time evolution of the Joint probability distribution function of the state 
space (2nx 1) column vector x and the system parameters. 

The use of the Liouville equation in mechanics and statistical mechanics is 
of long standing and goes back to Maxwell (see for example [1,2,3]). These 
references do not consider random system parameters, and average quantities 
under equilibrium conditions is of main Interest. While not of direct interest to 
us, it is possible to adapt these early methods to our needs. We derive the 
needed form of the Liouville equation following a procedure suggested by 
Kozin [4] for systems with random parameters and random or deterministic 
initial values. 

We are interested in the linear equations of motion in the form given by 

(1.3) with f= 0: 

x = Ax , (2.1) 

where x is the (2nxl) colwumn vector whose transpose x T has the form x T = 

jxj, . . . , x n ;x lt . . . , x n | and A is the (2nx2n) matrix given by the first of 

(1.4) The vector x is the state space form for representing the system response; 
the components of x will be denoted by x k (t), k=l, . . . ,2n. The random 
variables in A are denoted by X x , . . . ,X m . However, since the Liouville 
approach applies to general nonlinear as well as linear equations, we consider 
the general system 

x k = g k (Xj, . . . , x 2 jjj X lf . • . , X m ; t) , k=l, . . . , 2n . (2.2) 

Let p(x x , . . . , x 2n ;Xj, . . . ,X n ;t) be the Joint probability distribution of the 

random quantities (x lf . . . , x 2n ; X x X n ). We define the characteristic 

function <f> as 
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6 = E 


exp 1 


2n m ■k 

£ 0 k x k (t) + £<£jXj 


, i = v^T . 


k=i j=i > 

The differentiation of (2.3) with respect to time gives 


(2.3) 


d±_ 

<9t 


E 


2n ( 2n m \ 

i £^(<0 expi EMk^) + £<Ap^j 
k=l l k=l j=l i 


(2.4) 


The use of (2.2) in (2.4) yields 


86 2n 

lr = i 

k=i 


(2n m \ 

£0 k x k (t) + E^iEj 


(2.5) 


Since (2.3) is essentially the Fourier transform of the Joint density function 
p(x 1 ,...,x 2n ; X^-.Xjjjjt), the inverse Fourier transform of (2.5) produces 


— = -£ 

* U 


2n 5(gjP) 


cbc; 


( 2 . 6 ) 


The solution of (2.6) for p is given by a suitable function of the 
independent integrals of the Lagrangian system 


dt 

1 


- dp 


dx, 


dx 


2n 


{ d&i dg 2n 

I • * • T 


dx j 


dx 


2n 


®2n 


(2.7) 


Let u lf . . . , u 2n be 2n-independent integrals of (2.7). Then we know that the 
general solution of (2.6) is 


P (x x x 2n ; X x , . . . , X m ; t) = h (u lf . . . , u 2n ; X x , . . . ,X m ;t) (2.8) 

where h is an arbitrary function whose form is determined by the Initial 
conditions on x. 


In particular, consider the case where X x , . . . ,X m are explicit random 
mechanical parameters which are independent of the integrals. Then (2.8) can 
be written as 

p(x x » • ■ • » x 2n ; Xj, . . . , X m ; t) = h x (u x , . . . , u 2n ; t)h 2 (X x , . . . ,X m ) , (2.9) 

where h 2 is the Joint density function of the parameters, and h x is the 

conditional density of the randomness of the initial conditions. To illustrate the 
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form, consider as an example the simple linear system defined as 


( 2 . 10 ) 


x(t) + u> 2 x(t) = 0 ) 

x(t 0 ) = x 10 , x(t 0 ) = x 20 f 
In this case, we would have (x 10 ,x 20 ) random initial values, and u>, a random 


parameter. The Liouville equation is simply written as 


dp 


dp 2 dp 

= ~ X 2^— + WX 1^— 
u\j OX x uX 2 


where (2.7) is simply 


dt 


dx. 


dx c 


OJX 


i 


dp 

0 


In this case, we would easily find 


p(x 1 ,x 2 ;w;t) = hJu^x^Xao.t.Vw), u 2 (x 10 ,x 20 ,t,t o ;u;))h 2 (u;) 


^20 


= h x x 10 coso;(t- t Q ) H sina>(t-t 0 ) , 

* - - 


- u?x 1Q s\noj(t^. t 0 ) + x 20 cosa;(t- t D ) jh 2 (a;) . 

Upon utilizing the fact that the initial values can be defined as 


( 2 . 11 ) 


( 2 . 12 ) 


(2.13) 


x io = x 1 cosa;(t 0 - t) H sinw(t 0 -t) = u 1 (x 1 ,x 2 ,t 0 ,t,u;) 


U) 


(2.14) 


x 20 = - wX|Sinw(t 0 -t) + x 2 coso;(t 0 - t) = u 2 (x 1 ,x 2 ,t 0 ,t,o;) , 
the final probability form for the simple oscillator (2.10) becomes 

p(x 1 ,x 2 ;o;;t) = hj (u^Xj.x^to.tja;), u 2 (x 1 ,x 2 ,t 0 ,t;a;))h 2 (a') . (2.15) 

For the higher order system, the exact form of p in (2.8) would be obtained as 
in (2.15). 

We note that the initial conditions may be deterministic so that hj is a 
product of impulses at the origin and at unity 
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h x (u x , . . • , UgjjJ t) — 5( Uj) ...<5( u n ) <5( u n _j. j l)...6(u 2n ) i (2.16) 
where 6(.) is the delta function. Thus (2.15) would become 

p(x x , • • • , x 2n ; Xj, . ■ . , X m ; t| = 8( u j) ...6( u n ) S( u n _|_ x — 1 ) ...8( Uo n ) h 2 (Xj, • • ■ > -Xm) ■ 

(2.17) 


Let us illustrate this process by the simplification of the example (2.10)-(2.15). 
Consider again the undamped one degree of freedom linear oscillator. Let oj 2 
be a random parametric value. For t o =0, 


X 2 

u, =x 1 cosu;t- — slnwt , 

OJ 

u 2 = oj X ! sincut 4- x 2 coswt . 
Assume that oj 2 has a discrete distribution given by 


m 

ti 2 (w 2 ) = E p i <5(w 2 -cu 2 ) . 


For x x =0, x 2 =l at t= 0, (2.17) becomes 


(2.18) 


(2.19) 


p(x x ,x 2 ,u; 2 ; t) = «5(u x )<5(u 2 -l)E p i<5(u; 2 -u; 2 ) 


= 8 


1 

|x x cosu;t- — sinwtj 6 (u;x x sina>t + x 2 cosa;t- l)E p i^ w2 ~ <jJ ?) 

1 ( 2 . 20 ) 

Let us determine the mean of x x to Illustrate a possible use for (2.20); we have 


m 


E{x x } = f dx 1 f dx 2 f doj 2 |x x p(x x ,x 2 ,a; 2 ;tj 
Straight integration of (2.21) (see [7]) yields from (2.20) 


m sinuqt 

E{x x }= E P i 


( 2 . 21 ) 


( 2 . 22 ) 


Wi 


Other illustrations, including damping, are given in [4,5]. We are frequently 
concerned with the moments of x. Let us show how (2.6) can be employed to 
obtain them. 


To keep the details simple, consider the linear damped one degree of 


freedom system with equation of motion. 
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x x = x 2 (= g x ) 

k c , . 

x 2 = x x x 2 (= g 2 ) . 

m m 


(2.23) 


Equation (2.6) now takes the form 


dp 3(g x p) d(g 2 p) 

— — H 1 = 0 . 

dt dx x dx 2 


(2.24) 


Assume the m, k, c are independent of the initial vector, with h 2 (m,k,c) the 
probability density function of these parameters, and write 

p(xj,x 2 t;m,k,c) = h x (u 1 (x 1 ,x 2 ,t,m,k,c), u 2 (x 1 ,x 2 ,t,m,k,c)) h 2 (m,k,c) (2.25) 

In this case, we could simply rewrite the function h x as 

hi ( u i(xi,x 2 ,t,m,k,c)) = p 1 (x 1 ,x 2 ,t|m,k,c) (2.26) 

since the parameters (m,k,c) are conditional for h x . Upon inserting (2.26) into 
(2.25) and then into (2.24), we obtain 

dPi 5(g 1 p 1 ) 5(g 2 p x ) 


dt 


+ 


5x x 


+ 


dXr 


= 0 


(2.27) 


where h 2 (m,k,c) has been factored out. Let us evaluate the conditional 
expectations 


E{x 1 |m,k,c}= m 10 (t) = f Jx 1 p 1 dx 1 dx 2 

E{x 2 |m,k,c}= m 01 (t) = / / x 2 p 1 dx 1 dx 2 
Differentiation (partial) of these equations with respect to time produces 

<9p x 

“l,0 = I f X 1— dX l dX 2 


(2.28) 


(2.29) 


m 


o.i = J7 X : 


d Pi 
dt 


dx x dx 2 


But, from (2.23), 


5p x 

dt 


d(gjPi) ^(g 2 Pi) 

5x x 


dxr 


(2.30) 


The substitution of (2.30) into (2.29), the employment of (2.23) and simple 
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integration by parts of the resulting terms on the right-hand side of (2.29) 
finally yields 

m 10 = m o,i ’ 

(2.31) 

• k c 

m 0 j = - — m 10 - — m 0 j . 
m m 

The same procedure will produce the equations for the conditional moments 
E{xf |m,k,c}, E{x| |m,k,c}, etc. We note that for the first conditional moments 
we could have taken the conditional expectation of (2.23) to produce (2.31); 
however, this procedure only applies to the first moments. 

We integrate the moment equations (2.31) to obtain the conditional 
moments as a function of time. On multiplying these moments by h 2 (m,k,c) 
and integrating over m,k, and c, we obtain the moments of x x and x 2 . 

It is clear from the above discussion that the Liouville equation will provide 
the exact solution for the Joint probability density function p(x p . . . ,x 2n ; 
X lf ...,X m ;t) in the absence of external forces provided the integrals 
u lf . . . , u 2n can be obtained. Further, it provides a straightforward method for 
determining the moments of x from which means and variances of x can be 
obtained. 

The Liouville equation applies when there are no external forces. We are 
interested in the case when external forces are present, of course. Let us see 
what can be done along these lines. 

The Fokker-Planck equation is the natural extension of the Liouville 
equation (see [6,7]). We confine our attention to the case in which the 
external force vector f can be obtained by passing gaussian white noise through 
a stable linear damped system. We have as equations of motion, conditional on 
M= m, K= k, and C= c, 
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dx x = x 2 dt , 


dx 2 


k c X 3 

— x,dt- — x 2 dt 4 dt , 

m m m 


(2.32) 


dx 3 = - /?x 3 dt + dB , x 3 = 0 at t =0 

where we have employed the differential notation in this case, set f=x 3 , and 
were dB is the Brownian motion increment with 

E{dB}=0 , E{(dB) 2 } = a 2 dt . (2.33) 

The last equation of (2.32) represents the fact that the excitation is obtained by 

passing a gaussian white noise through a linear first order stable filter. We 

notice that for the Ito system (2.32) x T ={x 1 ,x 2 ,x 3 } is a vector Markoff process 

that generates a Fokker-Planck equation. 


It can be shown that in this case the Fokker-Planck equation for the 
conditional probability density function p x is 

( 


dP! 


dt 


d 

<9x x 


(x 2 pj - 


d 

dx 2 



a 2 d 2 Pi 
2 dxi 



(-/?x 3 pj) 


(2.34) 

We observe that all but the last term on the right are the same as would have 
occurred In the L.iouville equation in the absence of f. Let the conditional 
moments be 


m k 1 ,k 2 ,k s 




= 1 1 f x i l x 2 2 X 3 S PiUj-x^Xg) dx 1 dx 2 dx 3 


(2.35) 


Then, proceeding as in the development of (2.31), we find 
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m i,o,o — m o,i,o 


m 


0,1,0 


m 


m 


1 , 0,0 


m 


m o,i,o + I — l m 


m 


0 , 0,1 


(2.36) 


m 0 ,o,i — ^m 0 , 0 ,i 

On multiplying the solutions of (2.36) by h 2 (m,k,c) and integrating out the 
condition in these three conditional moments, we finally obtain the moments of 
x as a function of time. We obtain in analogous fashion the differential 
equations for the second conditional moments; we do not do this as the steps 
are of a mechanical nature and not of direct interest. The main point to notice 
is that differential equations for the conditional moments of x can be obtained 
when an external force is present in the equations of motion provided this force 
is produced by passing white noise through a suitable filter. 

It is important to point out that for any gaussian external excitation the 
solution vector is gaussian conditioned on the random parameters. Therefore, 
all conditional moments can be obtained but not as easily as above [8]. 

The Liouville equation enabled us to obtain, in a straightforward manner, 
the exact expression for the conditional probability density function p. 
Reference to (2.34) suggests that it will be much more difficult to obtain Pj 
from this equation and we shall not pursue this line of thought further. 
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III. Perturbation Methods 

The references [9-17] address the eigenvalue (natural frequency) and 
eigenvector (normal mode) problem In structural systems by perturbation 
methods. Before discussing methods or techniques involved, it is important to 
understand at the outset that the geometry of the structure, how its equations 
of motion are assembled, the final mathematical form of the equations of 
motion, and how randomness in parameters is introduced have a profound 
influence on the nature of the results obtained. 

A structure’s geometry can be in the form of a linear array (chain) of 
elements that may, for example, consist of simple harmonic oscillators strung 
together in a line, beam segments continuously connected at a sequence of 
supports In a line, etc. The geometry Is the simplest possible in such 
arrangements. Plate or shell type structures have a two-dimensional grid-like 
geometry and are next in order of complexity. Finally, we have the general 
case in which one or two-dimensional geometries are interconnected in a 
complex manner. 

The equations of motion depend on the coordinate choice, particularly when 
the fact that mass is always distributed is taken into account. Reference [18] 
discusses methods of making this choice and illustrates the substantial 
difference in response that can occur due to different choices. Reference [16] 
also discusses a component mode synthesis method for selecting coordinates 
and assembling the equations of motion. A coordinate transformation of the 
equations of motion is sometimes employed as In [10,19] and the altered form 
of the equations may be advantageous. 

Let us briefly present here a typical perturbation procedure. We consider 
the free motion of a conservative system governed by the equations 
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Iy + Ky = 0 (3.1) 

where y Is the (nx 1) column vector with transpose 

y T = {y v . • • ,y n } (3.2) 

Now the elements in the symmetric stiffness matrix K are determined by the 
bars, beams, columns, Joints, etc., making up the structure, and the 
uncertainties in the structure reside in these elements. Let there be m 
structural elements, and let the stiffness matrix of the 1 th structural element be 

Kj = (l+X i )K j , 1 = 1 m , (3.3) 

which produces the (nxn) random stiffness matrix 

m 

K = E K, = {K )k } . (3.4) 

1 

The random variables X 1( . . . ,X m are regarded as small perturbation terms 
describing the uncertainty present in the structural elements and we assume 

E {Xj} = 0 , Var{X i }=a i 2 , (3.5) 

Kj Is the mean stiffness matrix of the I th element, K=K T , i.e., K is symmetric 

in the Kj k , Kj k is the random stiffness element corresponding to yj and y k , and 

we assume masses of the elements do not change. We note also that we can 

write (3.4) as 

K = K + E XjKj , K = EKj (3.6) 

which gives also 

E {K}= K and = Kj , (3.7) 

where K. Is the stiffness matrix of the structure with each member taking its 
mean stiffness. 

Assume normal mode motion 

y = a cos(u>t + 4>) 

with a the (nx 1) column vector defined by 


(3.8) 
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Oi • t C* n } • (3.9) 

Then, substituting (3.8) into (3.1), we obtain 

(K- o> 2 I)a = 0 , (3.10) 

where again I is the (nxn) unit matrix. 

The squared natural frequencies w 2 are determined by the n roots of the 
equation 


det |K - w 2 I | = 0 , (3.11) 

revealing that the cu 2 and u> v are random variables since K contains random 

variables. Let the random mode corresponding to u> T be the (nxl) column 

vector a T . Then we can write 


(K- w 2 I)o: r = 0 , 

with the usual orthogonality relations 

a: r T Ia! s = 0 , a r T Ka s = 0 if s 7 ^ r 

a r T Io; r = 1 , a r T Ka: r = u / 2 , 
where ”T” denotes transpose, as before. 


(3.12) 


(3.13) 


We are now interested in expressing the random variables cu r and a T in 
terms of a power series in the random variables X ( . Consider, for example, the 
r 1 * natural frequency u r of the system expressed in the form 

m mm 

w r = w r + Yj + XI S ^ +••• > (3.14) 

i=i i=i j=l 

where 04 represents the r th natural frequency of the mean system, and the 
Xj.Xij,... are to be determined. Once we know the Xj.Xjj, . . . , we can obtain 
statistical properties of oj t or any other quantity of interest. Let us consider a 
general formulation of this problem, considering natural frequencies and 
normal modes. We follow the method suggested by Zarghame [14] which 
appears well suited to computation. 
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Differentiate (3.12) with respect to 


Kj - 2a; r 


du> r 

~dX, 


+ (K- w r 2 I) 


da T 

~dX: 


= 0 


(3.15) 


Next premultiply (3.15) by a r T obtaining 


oc T 


K: - 2oj 


d u T 

r ax" 


I \ a r = 0 , 


(3.16) 


since by the symmetry of K(K=K T ) and (3.12) 

a r T (K- w r 2 I) = 0 . 


Thus, with the third of (3.13) 


doJ T 

ax. 


2u> 


1 « r T • 


(3.17) 


M “ w r 

This is to be evaluated at X 1 =...=X n =0 ( i.e ., X= 0); we obtain 


8uj t 


ax. 


1 aJKa, 


2uk 


(3.18) 


where the underbarred quantities are to be evaluated for the system with mean 
stiffness. We note that (3.18) gives the sensitivity coefficients [20,21] of oj t 
with respect to the X t . The importance of the sensitivity coefficients resides in 
the fact that they reveal by their magnitudes those that are either sensitive 
or insensitive to uncertainity in members values. 


The a v , k=l, n span the coordinate space; hence, we may write 


a« r 


= s 04S . 


(3.19) 


ax s j 

We substitute (3.19) into (3.15): 

(& - 2uj r l^-l) a r + (K - ajfl] £ = 0 . (3.20) 

Now premultiply (3.20) by aj, obtaining 

a )[ [Ki - 2w r |^-l)a r + a k T (K- W 2 l) £ PiP a ) = 0 * 



- 16 - 


or, for k^r and with the use of (3.13) 


(^k- “?)0h k) = ~ a k Kj^r • 

Differentiating the next to last of (3.13) with respect to Xj gives 

da r 


a r T I 


dX i 


= 0 


which on premultiplying (3.19) by a r T I then demonstrates that /?W = o. Thus, 

K,a r 


/?(.*) = - 
f-'n 




, ky^Y , 


(3.21) 


and hence from (3.19) 


da 


r _ , ^k 

~ — 2-i 


a 


3X ^ 2 2 ” k ’ (3.22) 

oXj k W] f- 

where the prime on E* means that k does not take the value r. When 

evaluated at X= 0, we have 


da r 


dX, 


a i T Ki «. 

= s'4 — 

j <£r ~ UL f 


(3.23) 


where again the underbarred quantities are evaluated when members take on 
their mean stiffnesses. We note that (3.23) gives the sensitivity coefficients of 
the mode shapes with respect to Xj. Without showing the detailed deviation, we 
simply state that It can be shown that 


d 2 uj T 


dx.dx j , 0 


1_ 

2 


E (^k 2 - ^ 2 ) - 2 


dui j ) ( du>i 


dXiU^i'o 


(3.24) 


and 


where 


d 2 a T 


dX { dX j 


O 


E £r (k) • 


(3.25) 
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£$■ 


s* 2 - w k 



[ du { \ ) 

da T ) r 

( du x \ l 

' <9a r ’ 


T 

«k 

kMaxJJ 



. <9Xj , 

0 . 


&<!/ = - 


<9a r T 




<9a r 


ax 


, k^r 

(3.20) 

(3.27) 


i ' o 


Summarizing our results up to this point we have for the random variable 


W r = ^r+E 


du x 


ax. 


0 Xi+ 2^ i axjaxj , 0 


a 2 cu r 


XjXj+..., 


(3.28) 


where the partial derivatives are supplied by (3.18) and (3.24). We also have 
for the random variable 


a r = QL r + E 


da ] I 

^-L X ‘ + 7^ 


a 2 a r 

ax.ax 


XjXj +..., 


(3.29) 


j ' o 


where (3.23) supplies the first partial derivative, and (3.20) and (3.27) supply 
the derivatives in the double sum. 

Let us now consider the statistics of cu r , etc. Consider Eq. (3.28) first. We 
have, on taking expectation, 

i ( a 2 cu, i 

ER> = * + 7 ES |^-) 0 B < x ' x i> + - • (3 - 30) 

Even if the X’s are independent E{oj x }^oj^, since the E{X 2 }t^O terms are still 
present. Now square (3.28) and take expectation 

d 2 u) x 


UK) = 

* i ] 


ax^Xj 


E{X t Xj} 


+ EE 


duj r 


8oj t 


axJoiaxj'o 


i 


+ -EEE 

* i j k 


+ EEEE 

i i k g 


du x 


dX 


\ ’ O 


EfXiXj} 
d 2 u> x 


dX j<9X k , 0 


<9 2 u; r | 


a 2 w r i 

dX l dX j J 

0 

dX k dXy J 


ElX.XpCj 


EfriX^Xfi +. 


(3.31) 


We can now approximate Var u> r ; it is defined as 
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Var w r 2 = E{u; r 2 } - [E{u; r }] 2 . (3.32) 

Thus, it is a straightforward task to approximate the first two moments of ui T . 

If we extend (3.28) to cubic, quartric, ... terms in the X j( then (3.30) and 
(3.31) would contain additional terms. How far we should continue this 
process will depend on the relative size of the terms containing E{XjXj}, 
E{X,XjX k }, etc. and what information we have that would enable us to evaluate 
these expectations. It is not usual that we can evaluate any more than E{XjXj}. 

We note that Zarghame’s method described above gives sensitivity 
coefficients for natural frequencies and corresponding normal modes plus series 
expansions for these quantites In terms of the random variables X x , . . . ,X m 
that define the uncertainty present In the stiffness matrix K. Moments of the 
quantities are easily obtained, but it is practically Impossible to obtain 
distributions for the natural frequencies and corresponding normal modes. For 
confidence interval location and size for a natural frequency, for example, we 
must approximate using 

E{w r } ± 3 yVar (3.33) 

as a rough indication of a 99% confidence interval. This interval gives us some 

idea of the spread in a natural frequency and it could be employed to make 

reasonably sure that no steady excitation frequencies were contained therein for 

all uj t . Alternatively, we might employ the signal to noise ratio 

(3.34) 

yfv arUp 

to obtain an idea of how important stiffness uncertainty is for natural 
frequency; if (3.34) is greater than 20 or 30 say, we would regard the location 
of cJ r as deterministic; on the other hand, if (3.34) is less than 5-10, it might be 
unwise to ignore this level of variability in the location of w r , depending, of 
course, on the consquences of such uncertainty. 
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IV. Mean Square Approximate Systems 

We consider, in this section, a technique for including disorder or parameter 
uncertainty that follows a different line than taken In previous sections. 
Specifically, mean square systems are employed. We begin by introducing this 
type of system [22,23]. 

Let us begin with a very simple example in which there Is no disorder and 
no damping. Let the coordinates q 1( . . . , q n . Then, with 

2T=m jk q j q k , 2V = k jk qjq k , 5W=f j (t)6 qj , (4.1) 

where summation is on multiple subscripts. Then, with mass coefficients 

included, (4.1) can be rewritten as 

m jk Q k + k jk q k = fj(t) . (4.2) 

Let, with f 0 j constant, 

fj(t) = f 0j cos(u;t + <f>) . (4.3) 

Then, the forced motion 

q k = u k cos(o>t + <j>) (4.4) 

satisfies 

(k jk - w 2 m jk ) u k = f OJ . (4.5) 

These equations state that given the f 0 j and c u, the u k are determined by the 

solution of this linear system of equations. Further, if w is the natural 

frequency o> r and the u k define the r th mode shape o; rk , then the f oj must 

vanish. Let us look at the natural frequency problem in an unorthodox 

manner. 

Suppose we pick an ui and a set of u k which may not be one of the natural 
frequencies and normal modes. Then the right of (4.5) will not be zero and we 
need force amplitudes €j to produce this motion: 
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(k jlc - w 2 m jk ) u k = e, . (4.6) 

The 6j are the amplitudes required to maintain the assumed motion; we regard 

the 6 j as the amplitudes of the constraint forces required to produce the 

motion. 


Consider next a summation of second order amplitudes 

I(n,w) = fj ef > 0 . (4.7) 

l 

For a fixed u>, this is a positive definitive quadratic function of the u's. We can 
use this equation to find the natural frequencies cu r and corresponding normal 
modes a rk . Assume the u's are normalized in some manner (for example, 
u n =l, or better mj k UjU k =l). For fixed w, we find the minimum of l(u,cu)>0. 
Notice that if u>=u ; r the u's that produce a minimum are the a rk and 
l(a: rk ,w r 2 )=0, since the £j=0, J=l, . . . , n, in this case. It follows that if for a 
specified oj we find the minimum of I(u,u> 2 ) and this minimum equals zero, 
then this u is a natural frequency and the u that produces this zero minimum 
are proportional to the corresponding normal mode. Let us consider another 
interesting aspect of this method. 


Consider a frequency window g(w) with the following properties: 
g(w) >0 , oj' < oj < oj" , 

u" w" 

f g(uj)dtu = 1 , J u> 2 g(u/)doj < oo . 

0 O' U f 

Replace (4.7) with 


(4.8) 


w" 

I(u,g) = / ^e 2 du . (4.9) 

Find the the u that makes (4.9) a minimum. The interesting feature of this 
method is that if there is a natural frequency of the system in the frequency 
interval (o/,o>"), the u In minl(u,g) determine the normal mode of this natural 
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frequency. Let these u be In component form {uf r \ , , . ,u[ r '}; then the 
corresponding natural frequency is determined by the Rayleigh quotient: 

k jk u/ r) u( r) 


v? = 


(4.10) 


m;kUj (l M r) 

where we assume we have found the r 01 normal mode and its natural 
frequency. It follows that if there is concern that an Interval (u/,a>") contains a 
natural frequency, we have the method for determining if this is the case 
without determining all natural frequencies of the system. References [22,23] 
give details on this matter we will not discuss in this report. 


The computational problem of finding the minimum of I(u,o> 2 ) is carried 
out using one of a number of computer codes based upon conjugate gradient 
techniques, and, hence, is not a problem. 


So far, there has been no disorder in our system; i.e., the parameters mj k 
and kj k have been assumed to take definite values. Let us assume at this point 
that mass and stiffness contain random variables. We define, in this case, 


I(u,g) = E 


/ g(w) E e j 2d ^. - 

w' 1 


where, in vector-matrix form 


E ej* = u T (K-iiL 2 M) T (K-^ 2 M)u . 
i 

n 

Since E only operates on E e j 2 * n (4.11), we have 

i 


(4.11) 


(4.12) 


E | £e 2 J = E {u T (K-^ 2 M) T (K-££ 2 M)uj , (4.13) 

and w is a fixed parameter in (4.13). We assume the u are parameters to be 
determined. Thus, (4.13) takes the form 
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E | E e fJ = u T E {(K-w 2 M) t (K-w 2 M)} VI . (4.14) 

We note that (K- &, 2 M) T =K- 0 £ 2 M because of the symmetry assumed for K and 
M. In all events, means and second moments of K and M are all the 
Information needed to determine the expectation in (4.14). 

We than proceed as In the deterministic case, since I(u,g) has a 
deterministic form. 


To relate (4.7) to (4.13), all we have to do Is assume 

g(w) = <5(w-&.) , (4.15) 

where < 5 (.) is the delta function. The substitution of (4.15) into (4.11) yields 


I(u,w) = u t E j(K-w 2 M) T (K-w 2 M)j u . (4.16) 

This expression differs from (4.7) because of the assumed random parameters 
in K and M. If in (4.7), w is a natural frequency of the deterministic system, 
l(u,a;)=0. The l(u,a;)>0 in (4.16) because of the random parameters. Use of 
this fact has been made in [32] to obtain an estimate of the variance of natural 
frequency w r ; the formula is 


Var u T = 


I(u r ,w r 2 ) 


(4.17) 


4w; 


where cu r Is the r th natural frequency and u r is the corresponding normal mode 
for the system with mean parameter values. Monte Carlo simulation [24] 
reveals that (4.17) can be conservative and a correction is suggested. Equation 
(4.17) is easy to use since a minimum for I is not required. Further, (4.17) 
provides a much simpler method for estimating the variance of u> T than given in 
Section III. However, mean square approximate systems do not provide any 
information on E{cu r } or on variability in mode shape. Let us next consider 


how these systems apply to estimating frequency response with parameter 
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uncertainty present. 

We take the equations of motion in the form 

Mq + Cq + Kq = f . (4.18) 

The frequency response and Z (cv) are defined as follows. For the 

external force, 

f = <5 jr e iwt , (4.19) 

with r fixed at <5j r =0 for J^r, <5 rr =l, the component form of (4.18) is 

Qj = Z" 1 (o;)e iwt (4.20) 

which is exact. 

Suppose we try to approximate (4.20) with 

^ = /? jr e lwt , (4.21) 

where the ft jr are not known in advance. The equations of motion now are not 

satisfied and we must introduce constraint forces e i to bring about their 

satisfaction as 


(K jk - u, 2 M jk + iu,C jk ) /? kr - 5 jr = € jr . (4.22) 

From 


I(0.W) = E {f^ir} • (4.23) 

where asterisk denotes complex conjugate. This I is Just like (4.11) except the 
ft have replaced the u. We find the ft that make (4.23) a minimum, denote 
this ft by ft. Then, ft={ft lt . . . , ft n } is the mean square approximate to the 
Zjp^u;). It can be shown that if the system is deterministic (i.e., contains no 
random parameters) the ft are exactly the 

The tj r are complex, hence, the right of (4.23), when written out, Is 


E [ |"(Kj k - w 2 M jk - iwCj k )/? k i<5 jr ] (Kji w 2 Mji+iu>Cj|)/? lr - <5 jr ]} 
i=i L 


(4.24) 


- 24 - 


It follows that the minimum of (4.23) Is for the real and imaginary parts of /? kr . 
This added complication poses no additional computational problem [25,26]. 

The method also supplies an error criterion that makes it possible to Judge 
the accuracy of the /? kr . 

References [25,26] describe in some detail how the above technique can be 
applied to estimating the frequency response in a number of structures with 
specific attention being paid to numerical details of the computations. 
Reference [1] also describes how these techniques can be applied to the 
construction of a sequence of approximants for a complex system by starting 
from a highly constrained initial system and gradually relaxing the constraints. 
In these three references, extensive use Is made of the error criterion to 
determine when the estimated quantities (usually frequency response) are 
sufficiently accurate for the purpose In hand. A comment on what mean square 
approximate system provide is now in order. 

We observe, for example, that these systems enable us to estimate 
frequency response Zj^w) by means of /? kr (o;). The /? kr (u/) are deterministic 
numbers that take into account the means and variances of the statistical 
parameters of the structure. Thus, the /? kr (w) provide a deterministic estimate 
for Zj^a;). In the form given above and in [25,26,27], it is not possible to 
obtain statistical information concerning the However, it is possible to 

employ Monte Carlo methods to obtain estimates for the /? kr (w) given the 
parameters are sample values to obtain sample values for the /? kr (o;) from 
which statistical information can be obtained. 

The statistical energy approach (SEA) merits mention at this point since It 
also employs average energy concepts [28-32]. Basically, SEA estimates the 
average flow of energy from one part of a structure to another. For example, If 
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there is energy input into one part of a complex structure, this method provides 
an estimate of how this energy flows into another part of the structure. Thus, 
it is possible to estimate average vibrational energy present in any part of the 
structure. Information of this type is frequently the only type of information it 
is possible to obtain about the response in an extremely complex structure 
containing a large number of undamped natural frequencies in 1 Hertz. Insofar 
as we know, nothing has yet been done to include the influence of statistical 
parameters; however, the work given in [31] suggests it might be possible to do 


this. 
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V. Statistical Linearization of Nonlinear Systems 

1. Systems with Known Nonlinearity 

In a number of real systems the non-linearity property may be known. The 
statistical properties of the response of such systems cannot, in general, be 
determined. However, the statistical properties of the response of its linearized 
version may result in an adequate approximation. To obtain the linearized 
form, the method of statistical linearization has often been applied [33]. 

The basic idea is generated as follows. Consider the true system model 

y(t) + f(y(t)) = n(t) (5.1) 

where n(t) is the random excitation vector, which is assumed to be known as 

well. The idea is to approximate (5.1) as the linear form 

x(t) + Ax(t) = n(t) (5.2) 

The approach is to write (5.1) as, 

y(t) + Ay(t) + |f(y(t)) - Ay(t)j = n(t) (5.3) 

The linearization is obtained by choosing matrix A so that the ensemble 

average of the mean of the difference terms in (5.3) is minimized, i.e., 

mhi E | | |f(y( t) ) - Ay(t)|| 2 j (5.4) 

The solution is simply obtained as 

A = E (f(y)y') E {yy} -1 ( 5 - 5 ) 

where ( )' denotes transpose. 

Therefore, A is defined in terms of the y-statistics. But, the y-statistics are 
unknown! If they were known, we would not have to linearize the model. The 
traditional idea is to determine A via the linear X-system statistics as 

A = E {f(x)x'}E (xx'} _1 (5.6) 

The statistics in (5.6) can be exactly obtained in general, but when n is 
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Gaussian its determination defined by (5.2) becomes simpler. It should be 
pointed out that if f(y) is a vector polynomial then (5.5) can be determined 
exactly as statistical values, through stationary solutions of various forms [34], 

The following example for the Duffing oscillator with white noise excitation 
is very typical. 

Example I- Statistical linearization for the Duffing oscillator with white noise 
excitation. 

Consider the second order system 

y(t) + y(t) + y^t) = B(t) (5.7) 

The B process is the gaussian white noise, given through the classic Brownian 

process, which is Gaussian. The classic Ito form of (5.7) is written as 

dy x = y 2 dt 

■ , (5.8) 

dy 2 = - (y 2 +yi)dt + <*b 

From the conditions of B, we have 

E{dB} = 0 , E {(dB) 2 } = cr 2 dt (5.9) 

We wish to determine the linear form 

x(t) + x(t) + kx(t) = B(t) (5.10) 

which we can write as 

C dXj = x 2 dt 

1 dx 2 = - (Xa+kxjdt + dB t 5 - 11 ) 

Clearly, k is determined as 

k = E {x? • xjE {x 2 }- 1 = E { Xl 4 }E {x 2 }- 1 (5.12) 

However, since (x lt x 2 ) is gaussian, then we find 

k = 3E{x 2 } 

For the system (5.11) we can determine the stationary density as 


(5.13) 
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, ^ Vk 

p(x 1 ,x 2 ) = ex P 


7T <7 ' 


-T ( x l+ kx i ) 

a i 


This allows us to obtain the value of k and the x-moments as 


(5.14) 


k=1.22 a; E{x 1 x 2 } = 0 ; E{x 2 } = 0.408 a ; E{x|} = 0.500 a 2 (5.15) 

It is interesting to point out that for this particular example the moments of 

(y lt y 2 ) can also be obtained for the stationary case as 


E{ yi y 2 } = 0 ; E{y 2 } = 0.478 a ; E{yf} = 0.500 a 2 
We see that the error in the linear form is only in E{xf}, and is 


E{y I ^E{x 2 } 

E^ 2 } 

This is not an unacceptable error. 


0.070 

0.478 


or 14.4% 


(5.16) 


(5.17) 


2. Systems with Stochastic Parameters 

The question that must be considered here is, if the system has not only 
nonlinearities, but also possesses parameters that are random as well, in what 
way may we apply statistical linearization so that the random constant 
coefficients of the true system are reflected in the approximate system. We 
illustrate the procedure by considering again the Duffing equation 

y + cy + by 3 = B(t) (5.18) 

where here b is assumed to be a random constant. The linear form Is taken as 


x + cx + obx = B(t) 

where a, the linearization constant, is determined from 


which leads to 


min E j Hby 3 - ctby 1 1 2 
a ' 


a = E{b 2 y! 4 } E{b 2 y 2 } 1 

In the terms of conditional moments (5.18) can be expressed as 


(5.19) 


(5.20) 


(5.21) 
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a = E Jb 2 E{y x 4 |b}J E jb 2 E{y 2 |b}J 

Following the traditional path, we determine a from the linear x 
statistics, so that 


a 


= E jb 2 E {x(* |b}j E jb 2 E {x 2 |b}j 
= E|b 2 E{x 2 |b} 2 } E |b 2 {x 2 |b}} 


The conditional second moments are readily obtained as 

E{x,x 2 |b} = 0 ; E{x ? |b} = ; E{x| |b} = ^ 

which, when substituted into (5.23), leads to 

2 _ 3_ 1 _ 3_ a 2 

2 c E{b} 2 cb 

when the mean of b, E{b}=b. The x-moments can then be determined 


E{x 2 } = E |e{x 2 |b}J = 0.408(b/c) 1/2 a E{l/b} 

E{x 2 } = E |e{x| |b}J = 0.500 a 2 /c 

For illustration, assume a uniform distribution function for the 
constant b with mean b, 

f l/2eb , b(l-e)<b<b(l+0 
p ( b ) = l 0 , otherwise 


Thus 


E{l/b} = £n 
2eb 

where e=bandwidth parameters. Hence, 


( 1+0 

( 1-0 


E{x 2 }= 0.204(1 /be) 1/2 (<x/0# n 


(5.22) 

-system 


(5.23) 


(5.24) 


(5.25) 


(5.26) 

opening 


(5.27) 


(5.28) 


(5.29) 


The exact solution for this case is 
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E{y 2 } = 0.478(l/bc) 1/2 (<7/e) [( l+e) 1/2 - (l-e) 1/2 j 
and the error is 

E{y 2 }- E{x 2 } = Q.478 [(l+6) 1/2 -(l-£) 1/2 l - 0.204 gn(l+e/l 
E{y 2 } 0.478 [(l+e) l/2 -(l-e) 1/2 ] 

which for e=0 is 14.5% 


(5.30) 
— (5.31) 
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VI. Closure 

The requirement for precision pointing of large scale structures, consisting 
of trusses, radio reflectors and optical systems, has motivated the investigation 
of how parametric uncertainties affect the system response characteristics. 

Randomness of these large structures may arise from many sources, for 
example: manufacturing processes; space assembly by humans, variations in 
environmental conditions that may bear directly upon the material properties 
and resulting mechanical behavior. These uncertainties must be accounted for 
on some rational basis so that the quantities of concern for describing the 
response of the structure can be statistically characterized. Such quantities 
include: natural frequencies, normal modes, frequency response, etc. Specific 
random techniques that are available and discussed in this report are: 
Liouville’s equation, perturbation methods, mean-square approximate systems, 
and statistical linearization. 

These are the techniques that we shall consider in order to develop 
procedures for determining the response of space structures. It is conceivable, 
however, that advances in these techniques will have to be developed as our 


research effort unfolds. 
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